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^ Abstract 

-4— > Based on ab initio calculations and met adynamics simulations, we predict that 2H-M0S2, a 

g 

J ; layered insulator, will metallize under pressures in excess of 20-30 GPa. In the same pressure 

d 

^ range, simulations and enthalpy optimization predict a structural transition. Free mutual sliding 

of layers takes place at this transition, the original 2Hc stacking changing to a 2H(j stacking typical 

O of 2H-NbSe2, an event explaining for the first time previously mysterious X-ray diffraction and 

Raman spectroscopy data. Phonon and electron phonon calculations suggest that pristine M0S2, 

^ once metallized, will require ultrahigh pressures in order to develop superconductivity. 

00 

O PACS numbers: 61.50.Ks, 71.30.+h, 71.15.Pd, 74.62.Fj 

O 

o 

m 

> 

X 



1 



Transition metal dichalcogenides MX2, where M is a transition metal and X=S, Se and 
Te, are layered materials displaying a variety of electronic behavior from insulating to charge- 
density-wave to metallic to superconducting, with a rich scenario of phase transitions as a 
function of external parameters. Intercalation of electron-donating atomic species between 
the weakly bonded MX2 layers is facile, turning the insulators to metals which are interesting 
catalysts, and which also often superconduct at cryogenic temperatures. Perhaps the best 
known case in point is that of M0S2, a semiconductor with an indirect electronic band gap of 
1.29 eV pp, which in its pristine state is well known as a lubricant [2], and more recently as 
a potential photovoltaic [3] and single-layer transistor |1] material. M0S2 metallizes under 
intercalation of alkalis and alkaline earths [5], achieving superconducting temperatures up 
to 6 K. It has very recently also been shown to superconduct at even higher temperatures 
upon electric double layer (EDL) field doping [6l[7j. 

High pressure provides, besides doping, another general route to transform insulators to 
metals. Recently it was shown that bilayer sheet of M0S2 undergoes semiconductor-metal 
transition upon vertical compressive pressure [S]. The possibility that M0S2 might metallize 
under pressure is also suggested by a negative pressure coefficient of resistivity [9] indicating 
gap shrinking, dEc/dP < 0. Whether and how high pressure gap closing and metallization 
in bulk could be achieved is, however, an open question, also because X-ray diffraction [10] 
and Raman spectroscopy [11] data indicate a structural transition taking place between 20 
and 30 GPa to another phase of unknown structure and properties. What is the structural 
transition, whether high pressure metallization and possibly superconductivity will eventu- 
ally occur or not, and what could be the interplay among these structural and electronic 
phenomena in an initially insulating material, are all open questions of considerable interest, 
for which M0S2 can play a prototypical role. We alternated and combined well established 
first principles density functional theory (DFT) calculations and ab initio metadynamics 
(AIMtD) simulations [121 [13] to predict and clarify the simultaneous evolution of structural 
and electronic properties of pristine M0S2 under pressure. 

AIMtD simulation - where the supercell parameters act as collective variables - is a 
powerful computational tool to discover new and competing solid phases [HI [15]. For ab 
initio molecular dynamics and structural relaxations we employed VASP (Vienna ab initio 
simulation package) [TB] with standard scalar relativistic PAW pseudopotentials (no spin- 
orbit included) [T7] and cutoff of 340 eV. MD simulations were performed using a 72-atoms 
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2a/(3) X 2a/(3) supercell and 2x2x2 Monkhor st- Pack [T8] (MP) k-point sampling grid. For 
structural relaxations we used the 6-atoms unit cell and a 9 x 9 x 5 MP grid. The electronic 
structure calculations and structural optimizations were independently conducted using a 
variety of exchange-correlation functionals: LDA, GGA with PBE parametrization pjjj, hy- 
brid functionals B3LYP [20t |2T] and HSE06 [22]. Van der Waals interactions were included 
in the Grimme approximation [23] in the initial calculations at zero pressure, where they 
made an important contribution to the interplanar attraction, but were subsequently re- 
moved at high pressures. It turned out that at pressures above 5 GPa, in a regime where 
interlayer interactions are repulsive it is GGA with PBE potential and without van der 
Waals corrections that provides lattice parameters closest to experimental data. 

Within that approximation we thus carried out accurate total energy-based structural 
relaxations at increasing pressures, refining the properties of the phase(s) discovered with 
AIMtD dynamical simulations. That combined approach produced the results which we now 
describe. 




FIG. 1: Structures of 2Hc-MoS2 (a) and 2Ha-MoS2 (b). I. Brillouin zone for both structures [21] 
(c). 

Static relaxation of the initial structure 2Hc-MoS2 (space group PG^/mmc) of Figjl] (a) 
[32] upon increasing pressure produced no instability towards other phases up to 50 GPa. 
The unit cell dimensions shrank anisotropically (Fig|2]) as expected from a softer c-axis and 
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harder The indirect electronic band gap also shrank as expected, and eventually 

closed near 25 GPa with an electron-hole wavevector 5k (Figjsj) where K = (|,^,0) 

is a Brillouin Zone edge point (FigQc)) thus predicting band overlap metallization of 2Hc- 
M0S2 above 25 GPa. This result agrees with the very recent finding in Ref . [25] . Alternative 
DFT calculations and relaxations carried out with the HSE06 [22] and the B3LYP [20l 
[21] approximations yielded a metallization pressure of 35 and 55 GPa, respectively. Since 
these approaches, especially the latter, are known to underestimate metallicity, whereas 
GGA/PBE may overestimate it, we conclude that pressure induced metallization of 2Hc- 
M0S2 should be placed between 25 and 35 GPa. 

Extending in this pressure region our calculations to the 72-atoms 2a/(3) x 2a/(3) su- 
percell in place of the initial 6-atom unit cell, we further explored the possible onset of a 
periodic lattice distortion or charge-density- wave phase (more properly called an "excitonic 
insulator" [26]). We performed structural optimizations of samples with small random dis- 
placements of atoms and looked for the possible onset of the accompanying periodic lattice 
distortion, but failed to find any. In fact, the possible occurrence of an excitonic insulator 
driven superstructure in a narrow pressure interval near the insulator-metal transition, an 
occurrence which can be considered quite probable, cannot seriously be established within 
our methods, and more refined treatments of exchange will be called for in the future. 

We focused so far on the pressure evolution of the 2Hc-MoS2 initial structure. However, 
the local stability which we found upon relaxation is no guarantee for global stability at 
higher pressures. The AIMtD simulations, where the supercell vectors are biased away from 
the initial structure towards newer ones, can reveal the possible structures missed by straight 
relaxation. When applied to 2Hc-MoS2 in the 72-atom supercell at p = 40 GPa, it led to 
two successive and surprising sliding events shown in Fig|4j In the first event after metastep 
73 the two layers within the supercell shifted with respect to each other so that the Mo 
atom coordinates in the two layers became coincident (FigQa)). At the same time 

the supercell developed a tilt: the a and P angles moved away from right angle as can be 
seen in FigQc). In this intermediate structure the parameter a increased and the parameter 
c decreased from their original values (Figj4](b)). In the second event following metastep 
125 the a and /3 angles returned to the original 90 degrees while the parameters a and c 
further increased and decreased, respectively. The final structure was relaxed at p = 40 
GPa and T=0 and was identified as 2Ha - another 211 polytype with the same space group. 
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FIG. 2: Calculated lattice parameters a (up) and c (middle), and relative enthalpies (bottom) of 
2Ha- and 2Hc-MoS2 structures as a function of pressure. The experimental data taken from Table 
1 in Ref.[10] are in good agreement with calculated crossing of enthalpies at p ~ 20 GPa as shown 
in the bottom picture. 

but where the layer stacking is now AbA CbC in place of the initial AbA BaB stacking of 
2H^(27t [28]. In the new structure - typical of e.g., 2H-NbSe2 - all Mo atoms share the same 
(x, y) coordinates. Thus the pressure induced easy sliding of layers - a " superlubric" event, 
possibly related to MoS2's lubricant properties - transformed the 2Hc structure of M0S2 to 
2H„. 

Through accurate structural relaxations and electronic structure calculations we refined 
the cell parameters, the enthalpy and the electronic band structure of 2Ha-MoS2, with results 
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FIG. 3: Calculated band structures of 2Hc- (a) and 2Ha-MoS2 (b) at selected pressures. 



shown in Figs. [2| [3j The 2Hc and 2Ha enthalpies cross near 20 GPa, justifying the structural 
transition discovered by AIMtD.[33] The new 2Ha-MoS2 structure has at p = 20 GPa lattice 
parameters a = 3.080 A, c = 10.754 A with Mo atoms at Wyckoff position 2(b) and S 
atoms at 4(f) with z = 0.8947. As the band structures of Figjs] (b) show, the new phase 
2Ha-MoS2, also semiconducting at low pressures, undergoes a pressure-induced band overlap 
metallization now near 20 GPa. Structurally, 2Ha is less anisotropic than the initial 2Hc 
structure, the new interlayer spacing (reflected in the c axis length) smaller, the in-plane 
interatomic spacing (reflected by a, b lengths) larger. 

At this point comes an important contact with experiment. Few years ago Aksoy et 
al. [in] had reported an unspecified structural transformation taking place in 2Hc-MoS2 
between 20 and 30 GPa, without change of space group symmetry. When we overlap their 
measured lattice parameters a and c with our calculated ones in Figj2] the conclusion that 
their observed transformation is precisely the 2Hc to 213^ transition strongly suggests itself. 
To further nail down that conclusion we calculated X-ray powder patterns for 2Hc-MoS2 
and 2Ha-MoS2 and compared them in FigjSjwith the measured ones. Again, the agreement 
is quite good. The pressure induced onset of the 2Ha phase is heralded by the birth of a 
(104) reflection (indexed as "006" in Ref. |T0]) and a much stronger (102) reflection. The 
simultaneous demise of the 2Hc-MoS2 phase is signaled by the disappearance of (105) and 
the drop of (103) reflections. 

Raman spectroscopy data [TT] are also available for M0S2 in the pressure interval up to 
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FIG. 4: (a) Evolution of enthalpy in the metadynamics simulation during transformation from 
2Hc-MoS2 to 2Ha-MoS2 via an intermediate structure, (b) Evolution of unit cell parameters a, b, c. 
The parameters a, b are equal to the size of the 72-atoms supercell in respective directions divided 
by 2\/3. (c) Evolution of supercell angles a, (5. The angle 7 (not shown) fluctuates around 120° 
and does not undergo a significant change. 
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FIG. 5: Comparison of experimental [TO] and calculated X-ray diffraction patterns for 2Hc structure 
below 20 GPa and 2Ha structure above 20 GPa. 



31 GPa. At j9 = 19.1 GPa a new peak (called "d" band, Figs. 5, 6 in Ref.[TT]) appears, 
first as a shoulder on the high-frequency side of the original peak. As pressure grows its 
weight increases at the expense of the original peak and at 23.4 and 31 GPa both peaks 
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are clearly visible in the spectra. Using the Quantum Espresso code[29] we calculated the 
frequencies of Raman active and Ai^ phonon modes of 2Hc-MoS2 and 2Ha-MoS2 at the 
same experimental pressures of 19.1, 23.4 and 31 GPa.jM] Comparison with experimental 
data|Il]| in Fig|6] suggests that the "d" band coincides with the mode of the emergent 2Ha 
phase. On the other hand the difference of frequency of the higher Ai^ mode in both phases 
is very small which explains why this peak is not observed to split in the experiment. Apart 
from a small systematic difference of about 15 cm~^ between experimental and calculated 
peak positions the overall agreement is very good. This provides further strong support for 
the 2Hc — 7- 2Ha structural transition starting at 20 GPa with both phases coexisting up to 
the highest pressure of 31 GPa. 




FIG. 6: Comparison of experimental [llj and calculated peak positions of Raman spectra for 2Hc 
and 2Ha structures at different pressures. The vertical axis for theoretical data (right) has been 
shifted with respect to the axis for experimental data (left) by 15 cm~^ upwards to account for a 
small constant difference between the two datasets. 

We can now stretch beyond experimental observations, which stop just below 40 GPa, 
by extending simulations to higher pressures. |35] At 120 GPa, searching for newer potential 
structural transitions, we carried out another metadynamics simulation but the 2Ha struc- 
ture within the 72-atom supercell remained stable and did not show any transition. We 
further performed a structural relaxation of the 2Ha structure within the 6-atom unit cell 
at ultrahigh pressure of 200 GPa and again failed to find any new phase, indicating good 
local stability of 2Ha. 
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Having discovered that M0S2 turns metallic above 20-30 GPa, it is a relevant question 
whether it will superconduct. Inspection of the electronic density of states (DOS) at the 
Fermi level n{Ep) (Fig. [?]) shows that the metallicity of 2Ha-MoS2, marginal near 20 GPa, 
increases only very mildly with pressure. Poor metallicity is probably the reason why Raman 
spectra are still well defined even at 31 GPa. Using again the Quantum Espresso code|29j 
we calculated besides the DOS, also the dimensionless electron-phonon coupling strength A 
and the effective logarithmic average phonon frequency uiog at several pressures. As Table 
[l] shows, superconducting temperatures estimated by the Allen-Dynes formula |30J are not 
encouraging. Even for optimistic values for the Coulomb pseudopotential fi*, superconduc- 
tivity is predicted to appear only at ultrahigh pressures well beyond 100 GPa. 




E(eV) 

FIG. 7: 2Ha-MoS2 electronic density of states at selected pressures. 



TABLE I: Values of p [GPa], A, ujiog [K], niEp) [electrons/unit cell/eV], Tdfi* = 0.1) [K], Tdn* = 
0.05) [K] for 2Hc-MoS2 and 2H„-MoS2. 
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Our overall physical conclusions are therefore that pressure will cause M0S2 layers to 
slide, causing a structural transition between 2Hc and 2Ha polytypes, via a process akin 



to superlubric sliding. |3T] In both structures, band-overlap metallization will take place at 
relatively low pressures. After that, however, the broad cleft between valence and conduc- 
tion bands, both with large Mo d-hand character, does not shrink fast enough to produce 
a comparably large metallicity to that generated by, e.g., alkali doping, or by EDL field 
doping [HI |7j. High pressure M0S2 is predicted to remain semimetallic, and only prone to 
superconductivity well beyond 100 GPa. The contrast with the higher superconducting tem- 
peratures obtained by alkali and by EDL electron doping stems first of all from the smaller 
DOS achievable by pressure, but can also be related to a higher coupling of electrons relative 
to holes, and to a dimensionality effect in the EDL case. One aspect that remains to be 
explored is the possible occurrence, to be pursued with methods beyond those currently 
available, of excitonic insulator driven charge or spin density waves in a narrow pressure 
range close to the metallization pressure. If such a phase did exist, it is likely to be preceded 
and/or followed by a superconducting region. 

E.T. thanks T. Kagayama, Y. Iwasa, A. Shukla, and M. Calandra for discussions and 
exchange of information. Work in Trieste was partly sponsored by EU- Japan Project LEM- 
SUPER, and by Sinergia Contract CRSII2i36287/l. L.H. and R.M were supported by the 
Slovak Research and Development Agency under Contract No. APVV-0558-10 and by the 
project implementation 26220220004 within the Research & Development Operational Pro- 
gramme funded by the ERDF. Part of calculations were performed in the Computing Centre 
of the Slovak Academy of Sciences using the National Supercomputing Infrastructure sup- 
ported from Structural Funds of EU. L.H. acknowledges a predoc fellowship held at SISSA 
during part of this project. 



Electronic address: martonak@fmph.uniba.sk 



[1] L. Gmelin, Gmelin Handbook of Inorganic and Organometallic Chemsitry (Springer- Verlag, 

Berlin, 1995), vol. B 7,8,9, p. 16. 
[2] M. Dallavalle, N. Sndig, and F. Zerbetto, Langmuir 28, 7393 (2012). 

[3] K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz, Phys. Rev. Lett. 105, 136805 (2010). 
[4] B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti, and A. Kis, Nat. Nano. 6, 147 
(2011). 



11 



[5] R. B. Somoano, V. Hadek, and A. Rembaum, The Journal of Chemical Physics 58, 697 (1973). 
[6] K. Taniguchi, A. Matsumoto, H. Shimotani, and H. Takagi, Appl. Phys. Lett. 101, 042603 
(2012). 

[7] J. T. Ye, Y. J. Zhang, R. Akashi, M. S. Bahramy, R. Arita, and Y. Iwasa, Science 338, 1193 
(2012). 

[8] S. Bhattacharyya and A. K. Singh, Phys. Rev. B 86, 075454 (2012). 

[9] M. Dave, R. Vaidya, S. G. Patel, and A. R. Jani, Buh. Mater. Sci. 27, 213216 (2004). 
[10] R. Aksoy, Y. Ma, E. Selvi, M. C. Chyu, A. Ertas, and A. White, Journal of Physics and 

Chemistry of Solids 67, 1914 (2006). 
[11] T. Livneh and E. Sterer, Phys. Rev. B 81, 195209 (2010). 
[12] R. Martoiiak, A. Laio, and M. Parrineho, Phys. Rev. Lett. 90, 075503 (2003). 
[13] R. Martohak, D. Donadio, A. R. Oganov, and M. Parrinello, Nature Materials 5, 623 (2006). 
[14] R. Martohak, in Modern methods of crystal structure prediction, edited by A. R. Oganov 

(Wiley- VCH, Berlin, 2011). 
[15] R. Martohak, Eur. Phys. J. B 79, 241252 (2011). 
[16] G. Kresse and J. Furthmiiller, Phys. Rev. B 54, 11169 (1996). 
[17] G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999). 
[18] H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976). 
[19] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865-3868 (1996). 
[20] A. D. Becke, J. Chem. Phys. 98, 5648 (1993). 

[21] P. J. Stephens, F. J. Devlin, C. F. Chabalowski, and M. J. Frisch, The Journal of Physical 

Chemistry 98, 11623 (1994). 
[22] A. V. Krukau, O. A. Vydrov, A. F. Izmaylov, and G. E. Scuseria, The Journal of Chemical 

Physics 125, 224106 (pages 5) (2006). 
[23] S. Grimme, Journal of Computational Chemistry 27, 1787 (2006), ISSN 1096-987X. 



[24] A. Kokalj, Comp. Mater. Sci. 28, 155 (2003), URL |http : / /www . xcrysden . org 
[25] H. Guo, T. Yang, P. Tao, Y. Wang, and Z. Zhang (2012), URL |http : //arxiv . org/abs/ 
11208. 5941v2 '. 

[26] D. Jerome, T. M. Rice, and W. Kohn, Phys. Rev. 158, 462 (1967). 

[27] H. Katzke, P. Toledano, and W. Depmeier, Phys. Rev. B 69, 134111 (2004). 

[28] F. E. Wickman and D. K. Smith, The American Mineralogist 55 (1970). 

12 



[29] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. 

Chiarotti, M. Cococcioni, I. Dabo, et al., J.Phys.:Condens.Matter 21, 395502 (2009). 
[30] P. B. Allen and R. C. Dynes, Phys. Rev. B 12, 905 (1975). 

[31] A. Vanossi, N. Manini, M. Urbakh, S. Zapperi, and E. Tosatti, Rev. Mod. Phys., in press 



(2012), URL jhttp : //arxiv ■ org/abs/ 1 112. 3234[ 

[32] We adopt here the notation of Toledano et al.[57], while in earlier works such as, e.g., Ref.| 
this structure was denoted as 2H;,. 

[33] Zero-point motion and entropic contributions for both phases were estimated within the quasi- 
harmonic approximation at p = 10 GPa and T = 300 K. The difference between the two phases 
was smaller than 1 meV/atom and therefore these contributions were neglected in the figure. 

[34] We used the PBE PAW pseudopotentials Mo.pbe-spn-kjpaw.UPF and S.pbe-n-kjpaw.UPF. 

[35] Above 50 GPa we used the hard version of VASP PAW pseudopotentials. 



13 



